#include "field.hpp"
#include "meshTypeFunctions.hpp"
#include "bc_orthogonalDischarge_velocity.hpp"
#include "bc_orthogonalDischarge_waterdepth.hpp"
#include "util.hpp"

#include "testHelpers.hpp"
#include <catch2/catch_test_macros.hpp>
#include <catch2/catch_approx.hpp>

#include <cmath>

using namespace hms;
using Type = hms::boundaryCondition::Type;
using Catch::Approx;

TEST_CASE("Constant orthogonal discharge"){

	auto test = [](const auto& dummyMesh){
		using MeshType = remove_qualifiers<decltype(dummyMesh)>;

		auto mesh { uniMeshAs<MeshType>(3,4,1,1,true) };
		using BC_v = OrthogonalDischarge_velocity  <MeshType>;
		using BC_d = OrthogonalDischarge_waterdepth<MeshType>;

		const BoundaryPatch
			& left   { mesh.boundaryPatch("left") },
			& right  { mesh.boundaryPatch("right") },
			& bottom { mesh.boundaryPatch("bottom") },
			& top    { mesh.boundaryPatch("top") };

		Field
			v { mesh, 2},
			d { mesh, 1};
		
		scalar qBoundary;

		d.values() = 2;

		auto setPatchAndTest = [&](
			const BoundaryPatch& p,
			Array2s&& vDomain,
			Array2s&& vGhost,
			scalar dGhost
		){
			v.domainValues(p).colwise() = vDomain;
			v.setBoundaryCondition( BC_v{v, p, d, qBoundary} ).update();
			bool success {
				v.ghostValues(p).isApprox( vGhost.replicate( 1,p.size() ) )
			};
			if ( !success )
				std::cerr << fmt::format(
					"Patch \"{}\", target: ({}), actual:\n{}\n",
					p.name(), vGhost.transpose(), v.ghostValues(p)
				);
			REQUIRE( success );
			REQUIRE( d.ghostValues(p).isApproxToConstant(dGhost) );
		};
		auto setToThrow = [&](
			const BoundaryPatch& p, Array2s&& vDomain
		){
			v.domainValues(p).colwise() = vDomain;
			const auto& bc {
				v.setBoundaryCondition( BC_v{v, p, d, qBoundary} )
			};
			REQUIRE_THROWS( bc.update() );
		};
		SECTION("Divergent flow detection"){
			qBoundary = 1;
			setToThrow(left  , {-1, 0});
			setToThrow(right , { 1, 0});
			setToThrow(bottom, { 0,-1});
			setToThrow(top   , { 0, 1});
			qBoundary = -1;
			setToThrow(left  , { 1, 0});
			setToThrow(right , {-1, 0});
			setToThrow(bottom, { 0, 1});
			setToThrow(top   , { 0,-1});
		}
		scalar vres, dres;
		vres = 1.12948577306215;
		dres = 3.541434602718;
		SECTION("Subcritical inflow"){
			qBoundary = 4;
			setPatchAndTest(left  , { 1.8, 0.1}, { vres, 0.1}, dres );
			setPatchAndTest(right , {-1.8, 0.1}, {-vres, 0.1}, dres );
			setPatchAndTest(bottom, {-0.1, 1.8}, {-0.1, vres}, dres );
			setPatchAndTest(top   , { 0.1,-1.8}, { 0.1,-vres}, dres );
		}
		SECTION("Subcritical outflow"){
			qBoundary = -4;
			setPatchAndTest(left  , {-1.8,-0.1}, {-vres,-0.1}, dres );
			setPatchAndTest(right , { 1.8,-0.1}, { vres,-0.1}, dres );
			setPatchAndTest(bottom, { 0.1,-1.8}, { 0.1,-vres}, dres );
			setPatchAndTest(top   , {-0.1, 1.8}, {-0.1, vres}, dres );
		}
		scalar vin { std::sqrt(2*phys::g) };
		vres = 1.52084676907187;
		dres = 5.58899172017659;
		SECTION("Critical inflow"){
			qBoundary = 8.5;
			setPatchAndTest(left  , { vin, 0}, { vres, 0}, dres );
			setPatchAndTest(right , {-vin, 0}, {-vres, 0}, dres );
			setPatchAndTest(bottom, { 0, vin}, { 0, vres}, dres );
			setPatchAndTest(top   , { 0,-vin}, { 0,-vres}, dres );
		}
		SECTION("Critical outflow"){
			qBoundary = -8.5;
			setPatchAndTest(left  , {-vin, 0}, {-vres, 0}, dres );
			setPatchAndTest(right , { vin, 0}, { vres, 0}, dres );
			setPatchAndTest(bottom, { 0,-vin}, { 0,-vres}, dres );
			setPatchAndTest(top   , { 0, vin}, { 0, vres}, dres );
		}
		vin = 5;
		vres = 1.35605682530193;
		dres = 5.89945778873888;
		SECTION("Supercritical inflow"){
			/* Should this even be allowed?
			 * The solver is using the values from the domain to compute the
			 * state in the ghost cell, but for supercritical inflow conditions
			 * in the domain, no information can physically flow to the ghost
			 * cells */
			/* lower values for qBoundary result in even lower Froude Numbers
			 * in the ghost cells */
			qBoundary = 8;
			setPatchAndTest(left  , { vin,-0.2}, { vres,-0.2}, dres );
			setPatchAndTest(right , {-vin, 0.2}, {-vres, 0.2}, dres );
			setPatchAndTest(bottom, { 0.2, vin}, { 0.2, vres}, dres );
			setPatchAndTest(top   , {-0.2,-vin}, {-0.2,-vres}, dres );
		}
		SECTION("Supercritical outflow"){
			/* supercritical outflow is probably fine */
			qBoundary = -8;
			setPatchAndTest(left  , {-vin,-0.2}, {-vres,-0.2}, dres );
			setPatchAndTest(right , { vin, 0.2}, { vres, 0.2}, dres );
			setPatchAndTest(bottom, { 0.2,-vin}, { 0.2,-vres}, dres );
			setPatchAndTest(top   , {-0.2, vin}, {-0.2, vres}, dres );
		}
	};
	TEST_ALL_MESHTYPES(test)
}